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ABSTRACT 


The phenomenon of dynamic fracture is of interest in a wide variety of contexts from 
high speed manufacturing and defence applications to geophysics and seismology. The 
formulation and analytical solutions of dynamic fracture problems are very complex. 
Experimental and numerical methods like finite element methods are used extensively 
for solving dynamic fracture problems. In present work an attempt is made to simulate 
the dynamic crack propagation in DCB specimen by using the Cohesive Zone Model 
technique. 

The interface crack problem is modeled by using cohesive zone theory. Cohesive zone 
theory takes into account various inelastic processes in the process zone existing ahead 
of crack tip. The cohesive zone model characterized by quadratic traction separation 
law is used to model the cohesive zone. A two dimensional plain strain formulation is 
employed for the analysis. For discretization of the cohesive zone interface cohesive 
elements are used. 

The finite element model is first validated with hypothetical data for isotropic DCB 
specimen. Results give smooth variation of energy release rate curve. This method is 
also applied to the experimental data on glass fabric/epoxy composites. Results shows 
that the energy release rate gives fairly smooth variation. 
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Chapter 1 


Introduction 

1.1 General Introduction 

Composite materials are the choice as well as the need of the present world. Their 
popularity is due to their inherent properties such as high strength to weight ratio, 
and most important, controllable anisotropy i.e. different strength levels can be 
achieved in different directions according to the need. 

Mechanically fastened joints are critical parts in composite aerospace applications. 
The material combination in joined parts, along with stress concentration around the 
hole, makes a joint very complex element to design. Predicting strength of these 
parts is of great practical interest. Hence it is worthwhile to investigate into failure 
modes in mechanically fastened composite structures. Composite material with 
mechanical fasteners is commonly susceptible to the following modes of failure: 


1 . Bearing 

2. Tension 

3. Shear out 

4. Pull out 

5. Delamination 







Figure 1.1: Failure modes of composite material 


Figure 1.1 explains the first five modes of failures. These are common failure 
modes for isotropic material as well as the composites or orthotropic materials. 
Composite materials show an extra mode of failure called delamination. It is 
considered to be of prime importance and researchers are still pondering on this 
failure mode. 

While making holes by drilling process, fiber debonding and ply separation are 
likely to occur. Fiber debonding leads to a crack in a lamina which can be 
developed into an inter-facial crack. Ply separation further leads to the inter-facial 
or inter-laminar crack, commonly referred as delammation. Thus, delamination 
generally starts at the hole boundary. In practice delamination can start from free 
surface, loaded surface or from foreign material or due to matrix decohesion. 
Delamination has distinguished effects on stress concentration near hole 
boundaries. It normally leads to stress relaxation and hence alters the stress 
distribution near hole boundary. Once delamination is initiated or in fracture 







mechanics terms, nucleated, it propagates along the interface to result in complete 
ply separation. Due to ply separations each ply bears the load individually. 
Naturally, this result into low load bearing capacity, compared to the composite 
laminate without flaw. This individuality accelerates the failure process and may 
lead to catastrophic failure. Hence this thesis is devoted to investigate delamination 
initiation and propagation. 

Delamination can be viewed as inter-facial crack, and hence, fracture mechanics 
concepts can be applied to investigate delamination. For this purpose cohesive zone 
theory is used to model delamination, as it has specific advantages over 
conventional fracture mechanics methods such as LEFM, elasto-plastic, or J- 
integrad method. It should be noted that the cohesive zone models are no nlin ear in 
nature and hence nonlinearity is introduced into the problem. 

1.2 Literature Review 

Though considerable amount of work has been done to study the fracture 
phenomena through numerical methods, only rectangular beam (DCB) specimens 
were given special interest. Owen and Shantaram [26] started the case of finite 
element method to study dynamic crack growth. They studied DCB specimen and 
pipeline problem under transient loading. 

Nishioka and Atluri [27] investigated the crack propagation and arrest in high 
strength steel DCB specimen using moving singular dynamic finite element 
procedure. An edge crack in a rectangular DCB specimen was propagated by 
inserting a wedge and the results were compared with experimental data. In another 
work, Nishioka and Atluri [28] presented the results of generation and prediction 
studies of dynamic crack propagation in plane stress and plane strain cases. The 
studies were conducted by using FEM, taking into account stress singularity near 
the crack tip. Variation of dynamic stress intensity factor with time and variation of 
dynamic fracture toughness with velocity were studied and compared with available 
experimental results. 

For solving the crack problem cohesive zone model was used by many researchers. 
The cohesive zone idea was originally proposed by Barrenblatt [4] proposed 



Cohesive zone model as a possible alternative to the concepts of fracture mechanics 
in perfectly brittle materials. 

Dugdale [5] extended the concept of cohesive zone theory to perfectly plastic 
materials, by postulating the existence of a process zone at the crack tip. Dugdale 
applied the cohesive zone model to predict yielding of thin ideal elastic- plastic thin 
sheets containing slits. 

Needleman [6,7] was one of the first to propose polynomial and exponential based 
traction-separation laws for Cohesive zone model. Needleman analyzed a range of 
problems using Cohesive zone model, some of them were, particle debonding in 
metal matrix [6], tensile decohesion along an interface [7,8,9]. 

Xu and Needleman [10] further used the above models to study the void nucleation 
at interface between particle and matrix. 

Tvergaard [11], developed a similar model to that of Needleman. The formulations 
were identical for purely normal separation but Tvergaard coupled normal and 
tangential displacements to predict coexistence of three modes of crack. He 
modeled the decohesion of transversely staggered SiC whiskers in an aluminium 
matrix and subsequent void formation by normal and tangential separation. He also 
investigated fibre fracture and decohesion for a similar system with transversely 
aligned fibres. 

Tvergaard and Hutchinson [12] forwarded a bilinear type of cohesive zone model 
and showed relation between crack resistance curve and fracture process 
parameters. 

Li and Chandra [13] were the first to show the effect of form/shape of cohesive 
zone model in predicting the crack growth resistance and the size of plastic zone 
surrounding the crack tip, besides the effect of material parameters. They foimd that 
in addition to cohesive strength and cohesive energy, the form/shape of traction 
separation law of cohesive zone model plays a critical role in determining the crack 
growth resistance which clearly contradicts with the then common belief about the 
superfluousness of the shape. Chandra et al [14] showed the form dependency by 
evaluating the response in titanium matrix composites reinforced by silicon carbide 
fibres with two different forms (bilinear and exponential) of Cohesive zone model. 
El Sayed and Sridharan [15] used Cohesive zone model to analyze two cases viz., a 
DCS specimen and a compressed composite beam. They proposed to work in terms 
of stresses and strains of the material of cohesive layer instead of traction and 



displacement jumps. They employed a model proposed by Needleman [6] with a 
failure criterion based on the energy release rates in different modes. 


1.3 Overview of the work 

As mentioned in the previous sections, that delamination is the one of the important 
failure modes present in the mechanically fastened joints. Simulation of the 
delamination initiation and propagation is attempted in this work. Delamination is 
considered as an interface crack. Cohesive zone model is used to simulate this 
crack. 

For modeling the cohesive zone, special type of elements called cohesive elements 
are formulated. The cohesive elements are the four-node elements with zero area. A 
two dimensional plain strain formulation is adopted for the calculating stiffness and 
forces for the cohesive element. 

The composite material is assumed to bear a linear stress strain relationship. The 
composite laminae are modeled using the four noded isoparametric element. The 
quadratic cohesive zone model is used as the cohesive zone's governing traction 
displacement law. This cohesive zone’s constitutive equations introduce a material 
nonlinearity into the problem. 

The two dimensional plain strain formulation is adopted to simulate the crack 
propagation in DCB specimen. A hypothetical problem is considered for steel 
specimen in which dynamic force pulse is applied at the top of the specimen. 
Variations of energy release rate with various other parameters are studied. Along 
with these, experimental input data (Chowdhary [25]) are used to simulate crack 
propagation in graphite/epoxy DCB specimen. 

1.4 Thesis Organization 

The thesis report is presented in five chapters, as follows; 

Chapter2 discusses the various processes involved in the process zone and the 
cohesive zone model which is to be used during the analysis. 



Chapters describes the basics of two dimensional plain strain formulation for the 
composite material and the cohesive zone. It also describes the calculation of 
energy release rates. 

Chapter4 describes the results and discussion. 


Chapters describes the conclusion of the present work and the scope of future work. 



Chapter 2 


Theory of Cohesive Zone and Failure 
Criterion 


In this Chapter a brief discussion of various energies involved in cohesive zone and 
the model are presented. 

2.1 Introduction 

Conventional fracture mechanics, like linear elastic fracture mechanics, uses some 
idealized conditions which deviate from real situations. These theories characterize 
the crack by using single parameters such as stress intensity factors or energy 
release rates. These theories suffer from a disadvantage that near crack tip the 
stresses are singular, which is practically unrealistic. They assume the crack to be 
sharp, which may not always be the condition. The elastoplastic theory is somewhat 
generalized but it loses its usefulness when inelastic processes, such as plasticity, do 
occur near the crack tip, which are often observed in most of the materials. 

While conventional methods fail to explain the fracture problem fully, cohesive 
zone theory is formed to model various intricacies involved. In the cohesive zone 
theory, 

1 . The stresses near the crack tip can be bounded. 

2. Various inelastic, irreversible processes such as plasticity can be included. 



3. Various dissipating mechanism active in fracture process zone can be 
captured and the contributions from the mechanisms can be delineated with 
fine spatial resolution. 


2.2 Micro-mechanical aspects of the Cohesive Zone Model 

A typical case of a crack in a body is shown in Figure 2.1. OA denotes the crack 
and A is crack tip. The cohesive zone is ahead of crack tip and is of length A. The 
estimate of cohesive zone length is given as (Freund [17]), 
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Figure 2.1: A crack in a body 
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( 2 . 1 ) 


Where, amax is the maximum cohesive strength; Kiappi is stress intensity factor due 
to applied load. Considering a region in the vicinity of the cohesive zone, shown in 
Figure 2.2, representing a crack tip process zone. Process zone can be defined as the 
region within the separating surfaces where surface traction values are nonzero. 
This clearly means that, processes occurring within the process zone are accounted 
through traction displacement law. In the most of cohesive zone models, the 



traction separation relation for interfaces are such that, with increasing interface 
separation, the traction across the interface reaches a mayimnm then decreases and 
eventually vanishes, indicating complete decohesion. 



Figure 2.2: Crack tip with embedded process zone 



Figure 2.3: Typical traction separation curve of CZM and partition between intrinsic 

and extrinsic dissipation 



For example consider a model shown in Figure 2. 3. The significance of points A, B, 

C in Figure 2.2 and in Figure 2.3 is as follows: 

Point A corresponds to maximum separation 5sep which means that at point A the 
material is completely fractured and hence tractions are zero at that point. Thus, 
point A represents the actual crack tip, beyond which the material is fractured. 

Point B corresponds to the maximum traction. The separation at this point is 
denoted as 6max- Li and Chandra [13] called this point as a crack tip for representing 
the various quantities associated with the crack tip at this point. Point C corresponds 
to zero traction and zero separation. From this point theoretically the fracture 
surface starts to split. Hence, point C is usually referred as theoretical crack tip. Li 
and Chandra [13] referred the region stretching out from point B, i.e., peak cohesive 
stress to point C as forward region. The region stretching from peak stress, point B 
to point A is called as wake region 

2.3 Inelastic processes in cohesive zone 

During the fracture process, energy flows into the cohesive zone and region 
surrounding the cohesive zone. The rate, at which the inelastic processes develop 
within the cohesive zone, depends on the rates of dissipations caused by various 
mechanisms acting within the zone. Li and Chandra [13] has augmented that, since 
the energy flow is linked to micro-mechanics of fracture process, there should be a 
link between the active micro-mechanisms of a given material and that of traction- 
displacement curve. 

Before going further, it is necessary to investigate the various processes in the 
following regions: the forward region of traction-displacement curve (Figure 2.3), 
where inelastic damage processes are occurring due to current loading conditions, 
and the wake region, where elastic unloading is taking place. The micro-structural 
damage mechanisms in forward region are called extrinsic dissipation, which 
generally promote the growth of crack, while dissipative mechanisms in wake 
region are called intrinsic dissipation. Generally, consumed energy in the intrinsic 
dissipation is partly distributed in boundary material and rest goes into 
crack/cohesive zone. Intrinsic dissipation depends on the length of newly formed 



crack surfaces, while extrinsic dissipation is inherent property of material and 
doesn't depend upon geometric aspects. 

Figure 2.4 shows the summery of many processes in the wake region and the 
forward region [13]. It can be noted that all the processes under intrinsic dissipation 
heading, such as crack deflection, crack closure, fiiction, etc., are dependent on 
crack surface length. All the processes, such as micro-cracking, void growth and 
coalescence in ductile fracture, phase transformations, which are included under the 
extrinsic dissipation heading, are material dependent. 
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Figure 2.4; Energy dissipation mechanisms in the wake and forward regions [13]. 


2.4 Traction Separation law 


The basic concept of cohesive zone is introduced by Barenblatt [4] for predicting 
material behavior near crack tip region. The idea is extend by Dugdale [5] for a thin 
steel slits problem. After that various traction-separation laws are put forward. A 
summery of the various proposed CZM along with their application (as used by the 
proposers), is given in a paper by Chandra et al [14], in a tabular form. 

In the current analysis a cohesive zone model forwarded by Tvergaard [11] is used. 
Basically the model is postulated by Neddleman [8] and extended by Tvergaard 
[11]. The model used here is extrinsic type without rate dependency or history 
dependency. The cohesive zone model gives relation between tractions (Tn ; Ts ) 
and displacements ([Un] ; [Us]) in local directions, n, s as follows, 



where (Jmax is a measure of bond strength in normal direction and is a material 
property, X is a normalized quantity, i.e., 0 < >. < 1 which couples normal and 
tangential behavior and is given by 


A = 




kj 

V y 




y 


(2.4) 


where, 6„ , 5s are characteristic lengths in n, s directions, [Uj, ], [Us ] are the 
displacement jumps across the interface. In equations 2.3, 2.4, an, as are material 
parameters that help in defining the maximum traction in the respective directions. 
The maximum tractions in the n, s directions are Omax? cts^max- 



In pure normal separation ([Us] - 0), the maximum traction is CTmax- Total separation 
occurs at [un] = 5„ and work of separation per unit interface area is given by 

^ 1C ~ ^ max (2.5) 


Similar analogies can be stated for pure shear out separation. Thus, the work of 
separation per unit interface area is given by, 


G 


l!C 


= a 



( 2 . 6 ) 


In common practice, choice of Os, Omax is made and the parameters 6n, 6s are chosen 
such that the energy release rates are matched. But it is a common belief that energy 
release rates are the most important parameters and other parameters have little 
effect on the fracture process (El Sayed and Sridharan [15]). 



Figure 2.5: Normal Traction Separation Cohesive law 
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Figure 2.6: The tractions on the cohesive surfaces for normal separation 


It should be noted that in the above cohesive zone model interaction of the different 
modes is achieved by a parameter^. Displacement jump in any one direction results 
in change of traction values for all the three directions. The normal traction versus 
crack opening parameter for the normal separation is shown in Figure 2.5. In Figure 
2.6 an elemental surfaces under normal mode of separation are shown. Here the n, s 
directions indicate local directions. Note that, all the above calculations are in local 
co-ordinate system. 

The stiffiiess variation for a normal separation mode is shown in Figure 2.7. The 
stiffness terms will be derived later in the next Chapter. The stiffness decreases and 
becomes zero as X = 1/3 is approached. Later it decreases to achieve a minimum (or 
maximum negative value) at point X, = 2/3. After this point, stiffness value increases 
and becomes zero at X, = 1. The positive stiffness corresponds to the extrinsic or 
forward region, while the negative stiffness corresponds to the intrinsic or wake 
region. 




Fig 2.7: Stiffness versus Separation for normal separation 


2.5 Failure criterion 

In this Section a failure criterion for the cohesive zone model, described above, is 
discussed. As mentioned previously, the delamination is considered as the crack 
propagation problem. When the element fails, the crack is said to be propagated 
across the cohesive element. Two failure criteria are present in the literature. One of 
them is based on the energy consumed during the failure and its comparison with 
the energy release rates [15]. In the present analysis, a failure criteria based on the 
normalized parameter X is employed. % represents the normalized parameter, which 



combines displacement jumps across the interface, for the three modes of crack, 
into a single quantity. The failure criterion is stated as follows: when X reaches the 
value 1, the element is said to be failed. Thus, the cohesive element is capable of 
transferring the load when 0 < X < 1. But, when X reaches the unity value, the 
tractions and stiffness at the gauss points are set to be equal to zero. Once the 
element has failed, the tractions and stiffness values are set to zero for further 
analysis. 

Closure 

In this chapter, a brief discussion of cohesive zone model is presented and the 
various inelastic processes occurring in the cohesive zone are summarized. In the 
next chapter the implantation are presented. 



Chapter 3 


Finite Element Formulation 


3.1 Introduction 

In this chapter a finite element formulation of cohesive zone model as applied to 
composite material is described. As the formiilation for discretization of composite 
material can be found in many finite element text books, a very brief account of this 
is provided. HoAvever, a detailed description of finite element formulation of 
cohesive zone is provided. 

3.2 FEM formulation 

Equations that govern the dynamic response of a structural are derived by using 
virtual work principle involving work of external forces, internal, inertial and 
viscous forces for any small admissible motions. 

The finite element solution involves discretization of domain (Q) into suitable 
elements and approximating the field variable interior to the element in terms of 
nodal values using suitable shape function. 

FEM equation can be derived from the virtual equation for the body as 



( 3 . 1 ) 


\[5s^ {<j}dV + \{5qY + \[SqJkAqW 

0. o. a 

= \{sqY {P,}dv + Jter{/>,}</s 

« r, 

where 

{6q}, {Se}: Virtual displacements and corresponding strains 
{ 4' ^ Velocity and acceleration component 
{ a } : Stress component 
p : Material density 
k(i : Damping coefficient 
{Pb} ; Body forces 
{Ps} : Surface traction 
Fa : Surface traction boundary 

Displacements within element can be expressed in terms of its nodal 

values{Q}^®^ 

As 


kr=[A^r{2f (3-2) 

Strain-displacement relation can be written as 

w’=[Brter (3.3) 


Stress-strain relation or constitutive equation is given as 

[ar = = [Df (3-4) 


where 

[N] - Matrix of shape function 

[B] Strain displacement matrix 

[D] Matrix of material properties 



The present analysis is made for an undamped system and neglecting body forces. 
The stresses, strains and displacements are expressed in terms of nodal variables as 
given in equations 3.2, 3.3 and 3.4. This on substitution in Eqn (3.1) gives, 






e=l 


Av 


(3.5) 


where summation is over all elements, Ne- 

since {6q} is taken arbitrary, satisfying kinematic conditions, the equation can be 
written as. 


e=l \^^(.) j e=\ y \j<*) ^ 

(3.6) 


The final set of assembled equations can be written as, 

[KKe}+[JwM={f) (3-7) 

where 

[M] = 2 [M] , Global mass matrijt 

[K] = X [K-] > Global stiffness matrix 

{F} = X {F} > Global applied force vector 

and 

= ^Bf^^[Df^[Bf^\j\d4dTj , Elemental stiffeess matrix 

n<'> 

lp[Nf^^[NY\j\d§dTj , Elemental mass matrix (3.8) 

n<.) 

, Elemental force vector 



3.3 Finite Element formulation for cohesive element 



Figure 3.1: General domain with fixed external 
and time dependent internal boundary 

Considering a boundary value problem, as shown in Figure 3.1. A domain Q with 
boundary F has a time dependent boundaries (may be formed by the crack surface) 
Fi. Fi is viewed as an internal time dependent boundary. The following briefly 
gives the details of the FE modeling. 

For a no body force condition, the principle of conservation of momentum is given 
by 

o^ijj = 0. 

If ay is statically admissible stress field and 5ui is kinematiclly admissible 
displacement field, then principle of minimum potential energy leads to the 
equation, 

f^ij.j§UidQ = 0 (3.9) 

n 

where, 6ui are virtual displacements. It is to be noted that, in deriving above 
equation, a multiplication of displacement field 8ui with momentum conservation 
equation and symmetry of ay, is used. Also in this formulation, use of the total 
Lagrangian formulation with infinitesimal strain measure is made. 

Applying divergence theorem to equation 3.9, the left hand side will yield, 



[^ayjdu.do. = [^(Tyjdu^.da- +r2) 

= ) (3.10) 


Equation 3.10 along with equation 3.9 results in the equation, 

= |r,<5«,dr+ 


the various terms in the above equation show time dependence due to load nature 
and time dependence of internal boundary. Let the equilibrium state at time t + At 
is to be sought. Then, writing explicit time dependence of various terms in above 
equation, leads to 

|,a,(< + A/)*,(r + Ar>in= |:2-,(< + A/)5u,(/ + Ar)rfr+ ||,|r,(/ + A/K(' + A'V ['2 

(3.11) 

The incremental values for each parameters involved in the above equation are 
defined as follows. 


AT; =Ti(t + At)-T(t) 

AUj =Ui(t + At)-Ui(t) 
Aay=ai3(t + At)-aij(t) 
ASy= 8 y(t + At)-eij(t) 


Substituting above deflections of incremental forms in equation 3.1 1, the following 
equation is obtained: 



f[(a,j(t)+A(j„)68„(t)+0i,(t)846,]jn+|A(jj5A6ijdn 
n n 

= J[Ti (t + At)5Ui (t) + Ti (t)6AUi + ATi6AUi ]dr + 

r 


J[t, (t + 4t)6u, (t)+ T, (t)5Au, + AT,6Au,]dr, 

*" 2(0 


Considering changes in variational displacement field and traction to be small for 
a small time step At and neglecting their product term, the above equation can be 
restated as, 

Ja,j(t)5A8,jdQ + ]A0„8As,j<in = |(r(()+Ar, )dr+ |(r,(t)+Ar, Vr, 

n n r r^Cr) 


(3.12) 


Equation 3.12 represents the virtual work for a general body with initial and time 
dependent boundaries. 




Figure 3.2: Cohesive zone included between internal time dependent 
boundaries 

Having derived the virtual work expression in a general sense, focus is turned on 
the cohesive zone. Figure 3.2 shows a magnified view of the cohesive zone. In 
Figure 3.2, 5r2(t) represents the boundary of fracture zone on domain side with 

representing the normal to this boundary. dTo(t) represents the boundary of 



cohesive zone with representing the normal to this boundary. These boundaries 

are shown separately in Figure 3.2. It is clear from Figure 3.2 that the normals 
should be in opposite directions. Hence 


and 


5r:(()=ar.(0 


considering the above relations, the relation between the tractions acting on the two 
surfaces can be stated as, 


T,(n^Jt))=-T,(n,r.(t)) (3.13) 

Making use of equation 3.13, equation 3.14 can be restated for a virtual work of 
cohesive zone as, 


|AT,8Au,dr. + = - jT,{t)5Au,dr. - |o,(t)6A£,da (3.14) 

Tc a 

At this stage it is assumed that, in initial configuration, the cohesive zone has no 
volume and hence, completely described by the boundaries and the tractions applied 
on them. Hence at any arbitrary time t, the incremental form of equation 3.14 
reduces to 


fATibAUidr, = - fTi(t)5AUidr, (3.15) 

Th) r. 

With the aid of tangent stif&iesses, the incremental form of the tractions can be 
written as, 


tJ-,=k„4u,]+6T,'‘ 


(3.16) 



The second term in the equation 3.16 takes care of rate dependent traction 
separation law. If the rate dependency is neglected, as specified in the previous 
Chapter, the equation becomes, 

ATi=^A[Uj] (3.17) 

The kij terms in the above equation are obtained by differentiating equations 2.3, 3.4 
with respect to displacement jumps [uj] i.e. ky = 5Ti /5[uj ]. Hence for respective 
local n, s directions, ky are given as. 


"27^ 

^max 

(l-2A + A^)+- 

fkr 

2 

(-2 + 2A) 

UJ 




\ / 



(-2 + 2 ; i ) 


(3.18) 


(3.19) 


Note that above relations are valid for rate independent traction-displacement law. 
The variation of kn with separation parameter X, was given in a Figure 2.7. 
Substituting equation 3.16 into equation 3.14 the variational form is obtained as, 


jkyA[uj]^AUidr^=- (3.20) 

rdO C(0 


For simplicity, surface of cohesive zone is divided into two parts, viz. upper surface 
and lower line. So, the total boundary at any time is viewed as an addition of these 
two lines. We can write explicitly, what it means, in terms of equation as. 


r,W=r;(0 

where, letters u and 1, denote the quantities related to upper line and lower line 
respectively. 



Rewriting the incremental variational form of equation 3.20 in terms of upper line 
and lower line, the virtual work equation for a cohesive zone is obtained as, 


jk,j4[uj]5Au,dC + Jk,jA[uj]5Au,dr; 

CW r^(t) 

= - |T,(t)6Au,dC - Jt, 


rc"(t) 


ri(t) 


(t)6Au,dr; 


(3.21) 


3.3.1 Finite Element Discretization 

The cohesive surface is discretized using a four-node zero area 2-D rectangular 
elements. These elements are generally called as cohesive elements. Thus, the 
discretization of equation 3.21 in space is carried out by stating displacement field, 
on each surface, as a function of nodal displacements. 

Nodal displacements on upper and lower lines are indicated by u“ (u) and u“ (1) 
respectively, where u and 1 denotes the upper line and lower line respectively, 
superscript a denotes the node number and subscript i denotes the direction (n, s). 
Thus for an element shown in Fig 3.3, one has a e (1 ; 2) and i e (1; 2). 

Next, the assumptions for the upper and lower lines are stated clearly and separately 
as. 


The upper line: 


U; =N“u“(u) and AUj=N“u“(u) 
The lower line: 

Uj=N“u“(l) and Au; =N“u“(l) 


Where, N“ denote the interpolation functions or shape functions. With the help of 
above assumptions the displacement jumps, for the upper and lower lines, are given 
as. 



1 





[Uii]=0 


Interface element 2 


Interface element 1 


Figure 3.3: Crack opening displacement in mode I for a linear interpolated 

displacement field 

Upper Line: ' 

[m,. ] = M,. (u) - w, (/) = (u) - N“ti^ (/) 

a[m,.] = Aw,(w)- Am,.(/) = N“Aw“ {u)-N“Au“ (/) 

Lower Line: 

[«,]= »,(;)-u,(!<)= wx(0-iv“< W 

a[w, ] = Au, (/) - Aw, (m) = Aw“ {l)-N“Auf (u) 

Using the above assumptions, the equation 3.21 is rewritten as, 
jk jj [n ^ Au ^ (u) - N ^ Au (i)]n “ 5 Au “ (u)dr “ (t) + 

ro“(t) 

jkij [n^Au • (l)- N'^AUj (u)]n“5Au“ (u)(irj (t) 

ri(t) 

= - jTi(t)N“5Au“(u)(ir“(t)- |TXt)N“5Au“(l)drKt) (3.22) 



Assuming linear variation of displacement field on upper line as well as on lower 
line, the shape functions are defined as, 

N“ =(l/2)(H-^o) 

where, and take values +1 or -1 depending upon the local coordinates of 

nodes, a e (1; 2) corresponds to the node number. Same shape functions are 
assumed for upper as well as lower line. 

Factoring out the virtual displacements and separating out the equilibrium equations 
for upper and lower lines, two equations are obtained as, 

r“(0 r“(,) 

[au>{i) - (/jlcff-; (r) \T,{t)N - < (() 

ri(<) r '(0 


Above equations can be written in matrix form as. 


The above matrix equations can be rewritten as, 

-K,(u), K,Xu) 


AuJ(l)l ^ jFiO)] 
Au5(u)J |Fi(u)J 


(3.23) 


where, 

r;(<l 

K(!h- F,(')a'“<(') 

ri(<) 



To make better understanding of above equations, explicit expressions for Kij for 
directions n, s and for lower line are provided as follows. 



r;« 



T'At) 


The remaining terms, likeK“^Kp^,K“^ are taken as null. Kp^,K®^,Kf may be 
derived for the upper and lower lines separately. 

After substituting the expressions for Ky, the general form of equation 4.16 
becomes, 
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(3.24) 


where Kiower , Kypper and Flower , Fupper represent assembled stif&iess matrices and 
force vectors for the lower and upper line respectively. Auiower> Auupper are local 
incremental displacements vectors for lower and upper lines respectively. From the 
expressions in equation 3.24, it can be inferred that, basic differences for kn, ks for 
upper and lower lines come through difference between displacement jumps for the 
surfaces, as the shape functions are assumed to be same for the lines. But note that 
in the equations 3.18, 3.19 second power of displacement jumps are involved. 
Hence, essentially these terms are same for both the lines. Hence, 

K. upper ~ K. lower 

Scrutinizing the equations for the tractions in n, s directions (equations 2.3, 2.4) it 
can be stated that, tractions are exactly opposite to each other, as first power of 
displacement jumps are involved in their calculations. Hence it can be concluded 
that, 


F upper “F lower 



Considering the relations between the stif&ess matrices and force vectors for upper 
and lower lines, as inferred above, the equation 3.24 is simplified as. 
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— K 
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(3.25) 


3.3.2 Global assembly 

For an element shown in Figure 3.3, with linear variable displacement field i.e. a € 

(1,2), Ky (1) would be a 2 X 2 matrix. Thus two 2x2 matrices viz. K „ (1), Kj (1) are 
to be assembled to form a 4 x 4 matrix Kiower • 

Similarly, force vector of size 4x1, assembled from 2x1 sized force vectors Fn, 
Fs. 

Thus, from equation 3.25 cohesive element stif&iess matrix is of 8 x 8 size and 
cohesive element force matrix is of 8 x 1 size. 

Note that in the formulation, the local co-ordinate system is employed i.e. the 
displacement forms encountered are in local co-ordinate system. Before going for 
the calculations one has to transform the displacement vector in the global co- 
ordinate system, into the local one. Once the matrices and force vectors are formed 
using the local co-ordinate systems, then at the time of assembly they had to be 
transformed back in the global co-ordinate system. Mathematically this implies that, 

Kyz = T^K„sT Fyz = TF„s 

Keeping the transformations in mmd, v^e next move to global assembly of the 
cohesive element stiffiiess matrices and force vectors. It is interesting to know that, 
the global assembly of the matrices and vectors are simply carried out as, that for 
stiffness matrices and force vectors of any other usual finite element. 



3.4.1 Integration Algorithms 


Equations of equilibrium governing linear dynamic response of a system of finite 
element are given by Eqn (3.7). 

Eqn 3 .7 can be written as, 

Fi(t)+FE(t) = F(t) (3.26) 

where, F i ( t ) are inertia forces, Fi ( t ) = [M]{Q} , Fe ( t ) are elastic forces, 

F E ( t ) = [K]{Q}, all of them being time dependent. Thus, in dynamic analysis, the 
displacements {Q}, {8}and {F} are all functions of time. 

Mathematically, for small deformation A t can be assumed. Eqn (3.7) represent a 
system of linear differential equation of second order and the solution to the 
equation can be obtained by standard procedures of solving differential equations. 
The procedures for solving differential equations are divided into direct integration 
and mode superposition of which former is preferred in wave propagation problem. 
In this integration scheme, there are many different methods, which can be 
classified as ‘Explicit’ or ‘Implicit’ whose advantages and disadvantages are given 
by Bathe[29]. 

In the present work Newmark Integration scheme is used. 

3.4.2 Newmark’s Method of Integration 

Newmark’s method is an extension of ‘Linear acceleration method’, in which, 
linear variation of acceleration is assumed from time t to t + At. 

Velocity and acceleration approximated in terms of displacements are 

“{9)='{«}+l(l-'y)'{9)+'5“(9}l^ (3-27) 


(3.28) 



where a and 5 are the parameters chosen suitably to have accuracy and stability. 
Usually, a and 5 are taken as 0.25 and 0.5 respectively to get unconditional 
stability. As mentioned earlier, equilibrium at t + A t is considered along with these 
approximations. Solving equations 3.8 and 3.9, we get expressions for ' ^ *{q] and 
in terms of unknown displacements ' ^ '{ 9 } .These are then substituted in 
the following equation to solve for ' ^ . 

}='*"{-?'} ( 3 - 29 ) 


Entire method can be summarised as follows 

1 . Form global stiffness matrix [K], and mass matrix [M]. 

2 Initialise “{gf},® { 4 },° {^}. 

3 Select time step A/ and calculate integration constants 
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(3.30) 



2 la ; 


Ug = A?(l - S} ; a, = SAt 


For effective stiffness, [ii''] = [a:]+ a„[M] (3.3 1) 

4 For each time step : 

(i) Calculate effective load at t + At 

+ [wX«; {9} + a; {?}+ o; {«)) 0.32) 

(ii) Solve for displacement at time t + At by (3.33) 

(iii) Calculate acceleration and velocity at time t + At 

(3.34) 



3.5 Energy Release Rate and its Determination 


In fracture mechanics a crack can be characterized by four different parameters 

a) Energy release rate ( Gi Gu Giii ) is energy based and is also applied to 
brittle and less ductile material. 

b) Stress intensity factor ( Ki , Ku , Km ) is stress based. This parameter is also 
applied to brittle and less ductile material. 

c) J-Integral (J) which has been developed to deal with ductile material and can 
also be applied to brittle material as well. 

d) Crack tip opening displacement (CTOD) is displacement based and is 
developed for ductile materials. 

From computational point of view, stress intensity factor require special element so 
as to simulate stress singularity at the crack tip whereas J-Integral and energy 
release do not require any such modeling of stress singularity. In energy release rate 
approach, energy is foimd out in entire domain and in J-Integral; evaluation is done 
along a path far away from crack tip. Thus both these approaches adroitly avoid 
analysis close to the crack tip and don’t need any modeling of stress singularity. 

Griffith [30] outlined that in brittle fracture, the extension of a crack requires the 
formation of new surface, he reasoned with its associated surface energy, 
consequently, crack in a brittle solid should advance when the reduction of the total 
potential energy of the body during a small amount of crack advance equals the 
surface energy of the new surface thereby created. Most of the energy release, as 
the crack advances, comes from the parts of the specimen which are adjacent to the 
cracked surface. 

Two important parameters need to be considered 

1) How much energy is released when a crack advance. 

2) Minimum energy required for crack to advance in forming two mew 
surfaces. 

First parameter is measured in terms of energy release rate denoted by G. G is a 
function of crack size in general. Energy release rate is the amount of energy 



release per unit increase in area during crack growth. This energy is supplied by the 
elastic energy in the body and by the loading system. 

The energy requirement for a crack to grow per unit area extension is called as 
crack resistance and is usually denoted by symbol R. Crack resistance is sum of 
energy required to 1) form new surface and 2) cause an elastic deformation. 

Both the available energy release rate and crack resistance are important to study 
the possibility of crack becoming critical. Obviously, when the available energy 
release rate far exceeds the crack resistance, the crack starts to grow at high speed. 
In the present study , ‘energy release rate‘ is adopted as it is a more comprehensive 
concept. 

As the crack advances, 

1) Stiffness of the component decreases 

2) Strain energy in the component either increases or decreases. 

3) Work is done on the component if there is an application of load. 

4) Energy is being consumed to create two new surfaces. 

In the current analysis a cohesive zone model forwarded by Tvergaard [1 1] is used. 
Basically the model is postulated by Needleman [8] and extended by Tvergaard 
[11]. The model used here is extrinsic type Avithout rate dependency or history 
dependency. The cohesive zone model gives relation between tractions (Tn ) and 
displacements ([Un]) in local direction ‘n’ as follows. 



For a model where maximum cohesive strength Omax and critical displacement 5 are 
constants the amount of energy per unit area spent in the separation process , G, as 
described by this specific traction separation law is given by 



(3.37) 



In common practice, choice of (Jmax is made and the parameter 5n are chosen such 
that the critical energy release rates are matched with the known data. 

Energy release rate (G) is found by adding up the cohesive energy contributions 
from all N cohesive elements distributed in the finite element model Due to the use 
of variable Omax and 6 the cohesive energy spent in the N cohesive elements during 
crack growth has to be determined by integration 

G=Zfe¥(“.), W.w (3.38) 

<=1 

During the crack propagation analysis energy release rate at each time step is 
calculated by the equation 3.38. For the first cohesive element failure the cohesive 
parameters are chosen such that the critical energy release rate values are matched 
with the experimental data. And for the other cohesive elements the values are 
chosen such that energy release rate varies smoothly. 


3.6 Closure 

In this chapter, formulation of FEM equation is described along with integration 
algorithms for solving equilibrium equation governing linear dynamic response of a 
system of finite element. This chapter also describes some general concepts in 
energy and energy release rate. 



Chapter 4 

Results And Discussion 


4.1 Introduction 

In this Chapter results of the dynamic crack propagation analysis are presented 
using the method described in the previous Chapters. First the method is validated 
for dynamic case. Later results are also presented for the available experimental 
data for crack propagation. 

Various aspects of the results involve, 

• Crack propagation under dynamic input force pulse. 

• Crack propagation under time dependent input displacement (based on 
experimental data). 

4.2 Test Problem 

For the analysis purpose of investigating various aspects, a hypothetical problem of 
cantilever beam specimen is considered (Fig 4.1).In this test problem cohesive 
elements are inserted at the interface of the DCB specimen. The simulation is done 
for different combinations of cohesive zone materials parameters, viz., critical 
displacement jumps and maximum tractions for each element. 
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Fig 4.1 DCB Specimen 

Material and geometric properties are as follows 

• Modulus of Elasticity E=210 x 10® N/m^ 

• Poisson ratio n = 0.3 

• Length of specimen L = 0.03 m 

• Width of specimen W = 0.0025 m 

• Thickness of specimen B = 0.012 m 

• Crack length a = 0.015 m 

The DCB specimen is analyzed by using four-noded isoparametric elements. 
Details are as follows 


• Number of elements 120 

• Number of nodes 124 

• Size of element in X-direction 1 .0 mm 

I 

• Size of element in Y-direction 1 .25 mm 



The cohesive elements are analyzed by four-noded linear elements. Details are as 
follows 

• Number of cohesive elements 30 

• Number of nodes 3 1 

• Size of cohesive elements in X-direction 1 .0 mm 

For dynamic crack propagation, a short pulse is applied at the top of the specimen (Fig 
4.2). The pulse can be modeled by F( l-cos(cot) ) whose details are as follows: 

• Pulse duration p = 52 ps 

• Force amplitude Fo = 12N 

• Frequency of pulse f= 19.23 KHz 

Time step depends on the element size and wave velocity 



Fig 4.2 Input Force Pulse 


37 



Crack propagation starts after the force pulse of 52 ps is completely applied. Two 
different time steps are taken for the analysis 

1) Before crack propagation Ati = 4.92 x lO"^ sec 

2) After crack propagation initiates At 2 = a reduced value for higher 
accuracy. 

The meshing used for the DCB specimen is shown in Fig 4.3. 
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Fig 4.3 Meshing of DCB Specimen 

The Table 4.1 lists different cases of cohesive element material properties studied. 
For the first cohesive element it is assumed that the critical energy release rates are 
Gic= 512 J/m^ 

Giic = 904 J/m^ 

choice of Us , Omax is made and the parameters 5n , 5s are chosen such that the 
critical energy release rates are match well with the above mentioned values. 

And for the other cohesive elements these parameters are chosen in such a way that 
the energy release rate for propagation decreases gradually. No effect of the change 
in Os is observed , as expected, because this case represents pure mode I 
propagation. 


Elements 

^max (MP^) 

5n(mm) 

5s (mm) 

Os 

First 

130 

0.007 

0.0308 

0.4 

Second 

30 

0.005236 

0.0308 

0.4 

Third 

30 

0.004416 

0.0308 

0.4 

Fourth 

30 

0.0170442 

0.0308 












Fifth 

20 

0.004012 

0.000308 

0.4 

Sixth 

20 

0.00561 

0.000308 

0.4 

Seventh 

20 

0.00901 

0.000308 

0.4 

Eighth 

20 

0.01309 

0.000308 

0.4 


Table 4.1; Combinations of material properties of cohesive element 

. At each and every time step energy release rate (G) is found by adding up the 
cohesive energy contributions from all N cohesive elements distributed in the finite 
element model Due to the use of variable xjmax and 6 the cohesive energy spent in 
the N cohesive elements during crack growth has to be determined by integration 

G=i; i=i,w 

1=1 

where On is calculated by the traction separation law [Eq2.1] where the integration 
over each element gives the cohesive energy spent in each element. 

Above equation shows that energy release rate for propagation depends on the 
cohesive parameters viz. critical displacement jumps along the interface and 
maximum cohesive strength. So, by varying these parameters, fluctuations of 
energy release rate for propagation can be controlled. 

The variation of energy release rate (G) with time is shown in Fig. 4.4. 






Fig 4.4 Energy release rate for steel specimen 


At the crack initiation, energy is being consumed to create two new surfaces and due to 
that the energy release rate drops. As the crack advances, more and more new surfaces 
are formed and energy release rate shows a gradual decreasing nature. The present 
work follows the same trend for energy release rate curve. 

Variation of energy release rate is such that it starts increasing up to the crack initiation 
point and then shows a smooth variation with small fluctuations at the inter-element 
boundaries. 
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The velocity of crack tip is calculated by taking the ratio of element size in crack 
propagation direction and the total time to cross that element. And it is observed that 
the propagation time taken to cross the first cohesive element is more than the second 
cohesive element and so on. Fig 4.5 shows that the crack propagation initiates with a 
low velocity and then increases as the propagation takes place. Variation of crack 
velocity with crack extension is shown in Fig 4.5. 
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Fig 4.5 Crack velocity for steel specimen 
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4.3 Crack Propagation in Glass Fabric/Epoxy DCB Specimen 


The present method is applied to determine the dynamic Gi for the experimental data as 
reported by Chowdhary [25] . 

From the mechanics point of view fibre composites are among the class of material 
called orthotropic materials. Generalised Hook’s Law for 2-D orthotropic laminate is 
given by Agarwal [31]. 

(Jy = QijSjj (4. 1) 

i,j=l,2....,6 

where cr. are the stress components, Q^J is the stiffhess matrix, Bj are engineering strain 
components. 

For 2-D orthotropic material, 
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Ei^ and are elastic modulus in longitudinal and transverse direction 
G; is the shear modulus associated with longitudinal and transverse axes 
and o-i-i are major and minor Poisson ratios. 

Relation between four of the five elastic constants is 

E.f 

Following data is used as input for finding energy release rate. 

• Displacement at cantilever end as a function of time 

• Material data and specimen geometry 


Material and geometric properties are as follows 


• Material type 

• Material density 
Elasticity constants 


• Poisson ratio 


Glass fabric/epoxy laminate 
p=1825Kg/m^ 

El=26x10^ Pa 
ET=6.6xlO®Pa 

Glt=3.5x10^ Pa 

vlt^0.21 

vtl=0.05331 


Two different experimental data is used. 
Details of geometric data is given in Table 4.2 


(4. 7) 



Table 4.2 Details of DCB specimen 


Expt No 



B(mm) 

a (mm) 

1 

76 

3.8 

25.1 

40.3 

2 

76 

3.8 

24.8 

40.0 


L - Length of the DCB specimen 
W - Width of specimen 
B - Thickness of specimen 
a - Crack length 

In the experiment 1, time history of displacement at the cantilever end is knovra. Fig 
4.6 gives displacement applied at the cantilever end as a function of time. And this data 
is used as a input for finding dynamic energy release rate. 



Fig 4.6 Displacement of the cantilever end (Experiment 1) 
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For determining the energy release rate following cohesive material parameters are 
taken for the experiments 1 [Table 4.3]. 

For the first cohesive element, choice of Os , Omax is made and the parameters 5n , 5s are 
chosen such that the energy release rates are matched with the experimental data 
(Chowdhary [25]). For the other cohesive elements these parameters are selected such 
that the energy release rate varies smoothly. 


Elements 

Umax (MPa) 

5n(mm) 

8s (mm) 

Os 

First 

100 

0.00175 

0.0077 

0.4 

Second 

10 

0.0513366 

0.002156 

0.4 

Third 

10 

0.0275366 

0.000616 

0.4 

Fourth 


0.0309366 

0.000616 

0.4 

H9HH 

10 

0.0411366 

0.000616 

0.4 

BBHHI 

10 

0.0105366 

0.000616 

0.4 

Seventh 

10 

0.0037366 

0.000616 

0.4 


Table 4.3: Combinations of material properties of cohesive element for 
experiment 1 

Energy release rate variation for experiment no 1 is shown in Fig 4.8. Similar results for 
force release model [25] is shown in Fig 4.7. 
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Fig 4.8 Energy release rate for experiment 1 (Present study) 


At the crack initiation, energy is being consumed to create two new surfaces and due to 
that the energy release rate drops. As the crack advances, more and more new surfaces 
are formed and energy release rate shows a gradual decreasing nature. The present 
work follows the same trend for energy release rate curve. 

In the present work, the variation of energy release rate is such that it starts increasing 
with the time, reaches the maximum value at the crack initiation and then it falls 
slightly at the first inter-element boundary of the cohesive element. And after that it 
shows a continuously decreasing pattern with the small fluctuations at the inter-element 
boundaries. In the force release model energy release rate shows the same nature up to 
the crack initiation point but after that at the first inter-element boundary it falls to a 
great extent and thereafter it shows large fluctuations at the other inter-element 
boundaries. 
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The initial velocity of crack tip is calculated for first cohesive element and the final 
velocity is taken from experiment 1 (Chowdhary [25]). It is assumed that the crack is 
moving with constant acceleration during propagation. Crack tip velocity varies such that 
crack propagation initiates with a low velocity and then increases as the propagation 
takes place. 



Fig 4.9 Crack tip velocity for Experiment 1 
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Fig 4.10 shows the variation of energy release rate with the crack tip velocity. As the 
crack propagation takes place energy release rate decreases with the crack tip velocity. 
At the low crack tip velocity energy release rate drops more and gradually it reduces at 
high crack tip velocity. 



Fig 4.10 Energy release rate v/s crack tip velocity 
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Fig 4.11 shows the variation of energy release rate with the crack extension. As the 
crack propagation takes place, energy release rate increases within each elements and 
falls at the inter-element boundaries. Decrease in energy release rate at the inter- 
element boundaries reduces as the crack extension takes place as shown in Fig 4. 1 1 . 



Fig 4.11 Energy release rate v/s crack extension 
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Displacement (mm) 


In the experiment 2, time history of displacement at the cantilever end is known. Fig 
4.12 gives displacement applied at the cantilever end as a function of time. And this 
data is used as a input for finding dynamic energy release rate. 



0 0.00002 0.00004 0.00006 0.00008 0.0001 0.00012 

Time (s) 


Fig 4.12 Displacement of the cantilever end (Experiment 2) 
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For determining the energy release rate following cohesive material parameters are 
taken for the experiments 2 [Table 4.4]. 

For the first cohesive element, Os , Omax and the parameters 6n , 5s are chosen such that 
the energy release rates match well with the experimental data (Chowdhary [25]). For 
the other cohesive elements these parameters are selected such that the energy release 
rate varies smoothly. 


Elements 

^max (MPa) 

6n(mm) 

8s (mm) 

as 

First 

80 

0.000233 

0.001026 

0.4 

Second 

30 

0.00027166 

0.000924 

0.4 

Third 

10 

0.0001428 

0.000924 

0.4 

Fourth 

10 

0.00164628 

0.000924 

0.4 

Fifth 

10 

0.00504628 

0.000924 

0.4 

Sixth 

10 

0.01524628 

0.000308 

0.4 

Seventh 

10 

0.02884628 

0.00009548 

0.4 


Table 4.4: Combinations of material properties of cohesive element for 
experiment 2 

Energy release rate variation for experiment no 2 is shown in Fig 4.14. Similar results for 
force release model [25] is shown in Fig 4.13. 
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Fig 4.14 Energy release rate for experiment 2 (Present study) 


At the crack initiation, energy is being consumed to create two new surfaces and 
due to that the energy release rate drops. As the crack advances more and more new 
surfaces are formed and energy release rate shows a gradual decreasing nature. The 
present work follows the same trend for energy release rate curve 

For experiment 2 energy release rate curve for present method is shown in Fig 4.14, 
and it shows that the energy release rate increases till the crack initiates and then it 
slightly falls and then smooth variation of curve is shown with slight fluctuations as 
compared to the force release model [Fig 4.13]. 
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The initial velocity of crack tip is calculated for first cohesive element and the final 
velocity is taken from experiment 2 (Chowdhary [25]). It is assumed that the crack 
is moving with constant acceleration during propagation. Crack tip velocity varies 
such that crack propagation initiates with some low velocity and then increases as 
the propagation takes place. 



Fig 4.15 Crack tip velocity for Experiment 2 
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Fig 4.16 shows the variation of energy release rate with the crack tip velocity. As the 
crack propagation takes place energy release rate decreases with the crack tip 
velocity. At the low crack tip velocity energy release rate drops more and gradually it 
reduces at high crack tip velocity. 



Fig 4.16 Energy release rate v/s crack tip velocity 
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Fig 4.17 shows the variation of energy release rate with the crack extension. As the 
crack propagation takes place, energy release rate increases within the cohesive 
elements and falls at the inter-element boundaries. Decrease in energy release rate at 
the inter-element boundaries reduces as the crack extension takes place as shown in 
Fig 4.17. 



Fig 4.17 Energy release rate v/s crack tip velocity 
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Closure 


In this chapter results for dynamic crack propagation analysis in DCB specimen is 
presented and discussed. The conclusions drawn from the above discussions are 
presented in the next chapter. 
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Chapter 5 


Conclusions and Future Scope 


In the present work the analysis of double cantilever beam is carried out. The 
delamination front is simulated for these cases. A parametric study of cohesive zone 
model for the DCB specimen is carried out. In this chapter the conclusions of the 
work are presented and future scope for the work is commented. 

6.1 Conclusions 

The conclusions are drawn as follows: 

1 . The delamination can be modeled and simulated as a crack nucleating and 
propagating along the interface of composite material laminae. 

2. The delamination for a DCB specimen of glass epoxy composite material is 
simulated and is governed by mixed mode of crack propagation, with 
dominant mode I. 

3. Model gives fairly smooth and more stable variation of energy release rate. 

4. Experimental results reveal that unlike force release model energy release 
rate doesn’t drop to a low value in very few time steps but falls gradually 
which is a more practical result. 



6.2 Scope for Future Work 


1. The crack propagation model may be further modified by modifying the 
cohesive properties to get more stable and smooth variation of energy release 
rate. 

2. Cohesive zone model can be extended for the mixed mode analysis of composite 
materials. 

3. This method can also be extended to delamination analysis for 3-Dimensional 
crack propagation problems. 
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